Annual forest cover variation in Madagascar at protected areas and national levels (2000-2021)

1 Introduction

When an assignment results more difficult than expected, it is useful to codify, document and share the solution found. This is the philosophy that guided the foundation by KfW and partners of the Mapme initiative, the main tool that we use here: it was developed to enable policy-makers, funders, researcher and consultant that design, monitor or evaluate projects and policies to share their solutions, best practices in using spatial data for this purpose. To reciprocate with this virtuous approach, we summarize here some challenges that we encountered and that were not trivial to solve, as well as the workarounds we found. We hope it will help 1) the evaluator to correctly interpret the indicators provided, 2) future analysts to improve these calculations or reproduce them in other countries, and 3) the developers of the mapme package including me) to further improve this wonderful tool.

Our assignment was to calculate the forest cover and annual forest cover loss of all terrestrial protected areas in Madagascar and also at national level. These calculations are intended to serve as an input for a KfW ex post evaluation of its support to conservation in this country.

The output data is available in Results section 4, but we recommend reading the documentation surrounding it for an appropriate interpretation of the data. This works is complementary to a broader and more in-depth material developed for research and training purpose and available online.

2 Data

2.1 Data on protected areas

2.1.1 WDPA database on protected areas

Madagascar protected areas referenced in the IUCN protected area database (WDPA). It includes the spatial boundaries and information about the protected areas (see the variable meaning and codification here). We display this information below.

Code

# Load the R libraries needed for the analysis ---------------------------------------
# Install latest version of mapme.biodiversity
remotes::install_github("fBedecarrats/mapme.biodiversity",
                        upgrade = "always")
# Load required packages
library(mapme.biodiversity) # To get spatial data and compute indicators
library(tidyverse) # To easily handle data
library(lubridate) # To handle spatial data
library(writexl) # To export the result in exel files
library(sf) # To manage spatial data
library(geodata) # To get national boundaries
library(wdpar) # To get IUCN data on protected areas (WDPA)
library(tmap) # To display maps
library(gt) # To make nice tables
library(cowplot) # for combined graphs

# Get and display the WPDA data -----------------------------------------------------

# Get national boundaries
contour_mada <- gadm(country = "Madagascar", resolution = 1, level = 0,
                     path = "data/GADM") %>%
  st_as_sf()
# Get WDPA data  
WDPA_Mada <- wdpa_fetch("Madagascar", wait = TRUE, # using the wdpar package
                        download_dir = "data/WDPA") # store a copy locally
# Display the interactive map
tmap_mode("view") # for interactive mode
tm_shape(contour_mada) +
  tm_borders() +
  tm_shape(WDPA_Mada) +  
  tm_polygons(col = "IUCN_CAT", alpha = 0.6, title = "IUCN category",
              id = "NAME",
              popup.vars = c("Type" = "DESIG",
                             "UICN category" = "IUCN_CAT",
                             "Declared area (ha)" = "REP_AREA",
                             "Year of the status" = "STATUS_YR")) +
  tmap_options(check.and.fix = TRUE)

2.1.2 Limitations of WDPA database

We reviewed this data and compared it with another database, created for the association Vahatra by Goodman and colleagues in the framework of monographs they established for 98 terrestrial protected areas in Madagascar (Goodman et al. 2018). From this examination and comparison, we note that WDPA database presents serious limitations for Madagascar.

First: missing data. We can see that the information on protected areas creation dates, status and managing organization is often missing in WDPA.

Code

WDPA_Mada %>%
  st_drop_geometry() %>%
  summarise(`Total number of protected areas` = n(),
            `Missing IUCN category` = sum(IUCN_CAT == "Not Reported"),
            `Missing creation year` = sum(STATUS_YR == 0),
            `Missing managing organization` = sum(MANG_AUTH == "Not Reported")) %>%
  pivot_longer(cols = everything(),
               names_to = " ",
               values_to = "Nombre d'aires") %>%
  gt() %>%
  tab_header("Missing values in WDPA data for Madagascar") %>%
  tab_source_note("Source : WDPA (November 2022)")

Missing values in WDPA data for Madagascar
	Nombre d'aires
Total number of protected areas	171
Missing IUCN category	90
Missing creation year	51
Missing managing organization	163
Source : WDPA (November 2022)

Second: inaccurate data. An ongoing revision process by de Montalembert reveals the information on status and type in WDPA is also often inaccurate in WDPA.

Third: inaccurate geometries. The spatial perimeters of protected areas are sometimes inaccurate, as demonstrates the comparison below with the (apparently more precise) Vahatra data:

Code

if (!dir.exists("data/AP_Vahatra_shp")) {
  # Get data from Vahatra (copied on github to exemplify an issue)
  download.file("https://github.com/mapme-initiative/mapme.biodiversity/files/9746104/AP_Vahatra_shp.zip",destfile = "data/AP_Vahatra_shp.zip")
  unzip("data/AP_Vahatra_shp.zip",exdir = "data/AP_Vahatra_shp")
}

PAs_Vahatra <- read_sf("data/AP_Vahatra_shp/AP_Vahatra.shp") %>%
  rename(cat_iucn = cat__cn, creation = creatin, date_creation = dt_crtn,
         hectares = hectars, full_name = full_nm) %>%
  st_make_valid() %>% # some geometries are invalid
  mutate(an_creation = year(date_creation))
# There is no CRS in the data. We assume WGS 84 is correct
st_crs(PAs_Vahatra)<-"EPSG:4326"

# On harmonise les noms qui sont parfois notés différemment entre les sources
PAs_Vahatra <- PAs_Vahatra %>%
  mutate(nom_wdpa = case_when(
    nom == "Corridor Forestier Bongolava" ~ "Corridor forestier Bongolava",
    nom == "Ranobe PK32" ~ "Ranobe PK 32",
    str_detect(nom, "Ambositra-Vondrozo") ~ "Corridor Forestier Ambositra Vondrozo",
    nom == "Réserve deTampolo" ~ "Réserve de Tampolo",
    nom == "Bombetoka Beloboka" ~ "Bombetoka Belemboka",
    nom == "Ampananganandehibe-Behasina" ~ "Ampanganandehibe-Behasina",
    nom == "Forêt Sacrée Alandraza Analavelo" ~ "Analavelona", # vérfié sur carte : les mêmes
    nom == "Réserve speciale Pointe à Larrée" ~ "Réserve spéciale Pointe à Larrée", 
    nom == "Vohidava-Betsimalaho" ~ "Vohidava Betsimalao", 
    nom == "Anjanaharibe Sud" ~ "Anjanaharibe_sud",
    nom == "Iles Radama/Sahamalaza" ~ "Sahamalaza Iles Radama",
    nom == "Kalambatritra" ~ "Kalambatrika",
    nom == "Mananara-Nord" ~ "Mananara Nord",
    nom == "Kirindy - Mitea" ~ "Kirindy Mite",
    nom == "Midongy du Sud" ~ "Befotaka Midongy", # Vérifié sur la carte
    nom == "Montagne d'Ambre/Forêt d'Ambre" ~ "Montagne d'Ambre",
    nom == "Tsimanampesotsa" ~ "Tsimanampesotse",
    nom == "Pic d'Ivohibe" ~ "Ivohibe",
    nom == "Forêt Naturelle de Petriky" ~ "Forêt Naturel de Petriky",
    nom == "Tsingy de Namoroka" ~ "Namoroka",
    nom == "Réserve de Ressources Naturelle Mahimborondro" ~ "Mahimborondro",
    str_detect(nom, "Complexe Tsimembo Manambolomaty") ~ "Complexe Tsimembo Manambolomaty",
    nom == "Mandrozo" ~ "Zone Humide de Mandrozo",
    nom == "Paysage Harmonieux Protégés Bemanevika" ~ "Complexe des Zones Humides de Bemanevika",
    nom == "Nord Ifotaky" ~ "INord fotaky",
    TRUE ~ nom)) %>%
  arrange(nom_wdpa) %>%
  mutate(rownum = row_number())

# On ne garde que les aires de WDPA qui apparaissent dans Vahatra
WDPA_commun <- WDPA_Mada %>%
  filter(NAME %in% PAs_Vahatra$nom_wdpa) %>%
  filter(!(NAME == "Analalava" & IUCN_CAT == "Not Reported")) %>%
  filter(!(NAME == "Site Bioculturel d'Antrema" & IUCN_CAT == "Not Reported")) %>%
  filter(DESIG != "UNESCO-MAB Biosphere Reserve") %>%
  arrange(NAME)  %>%
  mutate(rownum = row_number())
       
# Cette fonction calcule la part d'un polygone incluse dans un 
# autre polygone et retourne un ratio entre 0 et 1
ratio_inclus <- function(x, y) {
  inclus <- st_intersection(x, y)
  ratio <- st_area(inclus) / st_area(x)
  return(ratio)
}

# On calcule la part des polygones Vahatra incluse dans les polgones WDPA 
V_in_W <- map2_dbl(WDPA_commun$geometry, PAs_Vahatra$geometry, ratio_inclus)
# Puis l'inverse
W_in_V <- map2_dbl(PAs_Vahatra$geometry, WDPA_commun$geometry, ratio_inclus)
# On fait un facteur des deux
recoupement_mutuel <- V_in_W * W_in_V
# Qu'on ramène dans les jeux de données d'origine
WDPA_commun2 <- bind_cols(WDPA_commun, V_in_W = V_in_W, W_in_V = W_in_V,
                         recoupement_mutuel = recoupement_mutuel) %>%
  arrange(recoupement_mutuel, rownum)
PAs_Vahatra2 <- bind_cols(PAs_Vahatra, V_in_W = V_in_W, W_in_V = W_in_V,
                        recoupement_mutuel = recoupement_mutuel) %>%
  arrange(recoupement_mutuel, rownum)

min_recoup <- WDPA_commun2 %>%
  filter(NAME != "Nosy Mangabe", NAME != "Sahamalaza Iles Radama") %>%
  filter(row_number() <= 10) %>%
  select(nom_wdpa = NAME, rownum) %>%
  mutate(source = "WDPA") %>%
  bind_rows(select(filter(PAs_Vahatra2, rownum %in% .$rownum), nom_wdpa, rownum)) %>%
  mutate(source = ifelse(is.na(source), "Vahatra", source)) %>%
  arrange(nom_wdpa)
tmap_mode("plot")
min_recoup %>%
  tm_shape() +
  tm_polygons() +
  tm_facets(by = c("nom_wdpa", "source")) +
  tm_layout(panel.label.size=2.5)

Fourth: overlapping geometries. Several protected areas referenced in WDPA overlap one with another, either completely (protected areas within protected areas) or partially (shared territories). In order to gauge the extent of these overlaps, we first calculate the sum of the surface areas of the protected areas recorded in the WDPA database, and then compared it to their total footprint, without duplication: 21% of WDPA area is duplicated (see here for more details). This means that indicators such as total area, forest cover or forest cover loss related to the WDPA polygons cannot be summed or averaged (or we would be counting twice or thrice the same areas).

2.1.3 Data from the Vahatra association

From this, we would recommend to use the Vahatra data instead of the WDPA data, here is a description. The data from Vahatra is available on a dedicated portal that makes it accessible under a creative commons licence (CC-BY).

Code

tmap_mode("view")
tm_shape(contour_mada) +
  tm_borders() +
  tm_shape(PAs_Vahatra) + 
  tm_polygons(col = "cat_iucn", alpha = 0.6, title = "Catégorie IUCN",
              id = "nom",
              popup.vars = c("Creation" = "creation",
                             "Creation year" = "an_creation",
                             "Area (ha)" = "hectares",
                             "Full name" = "full_name",
                             "Manager" = "gest_1")) +
  tmap_options(check.and.fix = TRUE)

Hovever, some protected areas included in WDPA are missing in the Vahatra database. See here for more details.

2.2 Forest cover data

There are several definitions of forest (For a review of prominent definitions, see Chazdon et al. 2016). We retain the FAO (2000) as a reference:

Land with tree crown cover (or equivalent stocking level) of more than 10 % and area of more than 0.5 ha. The trees should be able to reach a minimum height of 5 m at maturity in situ. May consist either of closed forest formations where trees of various storeys and undergrowth cover a high proportion of the ground; or open forest formations with a continuous vegetation cover in which tree crown cover exceeds 10 %. Young natural stands and all plantations established for forestry purposes which have yet to reach a crown density of 10 % or tree height of 5 m are included under forest, as are areas normally forming part of the forest area which are temporarily unstocked as a result of human intervention or natural causes but which are expected to revert to forest.

Note that, as reported below, we were not able to perform the computation with the 0.5 ha threshold and had to use a 1 ha threshold.

3 Methods

The computation is performed using the mapme.biodiversity package (Görgen and Bhandari 2022).

Note that we faced difficulties:

Updating the source data (solved): The data fetched by mapme.biodiversity includes the forest cover losses for years 2000 to 2020, but not for 2021. We copied the package and modified it to fetch the most recent data that includes forest cover loss from 2021. This modification will soon be submitted to the original mapme.biodiversity package.

Minimum forest patch size incompatible with FAO definition (unsolved): The mapme.biodiversity packages is programmed to define as forest areas of a minimum size, and this minimum size must be a full number in hectares (so at least 1 hectare). However, thhe FAO definition considers an area as forest with patches starting from 0.5 ha. We intented some workarounds that did not work, so we have to retain 1 ha as the minimal extent of forest patches for now. An issue has been filled to be added to the mapme development backlog.

Handling polygons that overlap satellite image tiles (solved): The source data fetched by mapme.biodiversity come in tiles (3 in the case of Madagascar). The package cannot compute indicators for polygons that overlap the tile borders, which is the case for the national territory of Madagascar, as well as for 7 protected areas (see image below).

Code

# The following code sequence launches very heavy data processing. We want to avoid
# re-launching it if not needed. So we start checking if the output of the data 
# processing is available locally, and only if not we launch it.
if (file.exists("data/outputs_GFC.rds")) {
  load("data/outputs_GFC.rds")
} else {
  load("data/outputs_GFC.rds")
  # Install latest version of mapme.biodiversity
  remotes::install_github("fBedecarrats/mapme.biodiversity",
                          upgrade = "always", force = TRUE)
  
  WDPA_poly  <- WDPA_Mada %>%
    filter(st_geometry_type(.) == "MULTIPOLYGON") %>% 
    st_cast("POLYGON")
  
  WDPA_poly <- init_portfolio(WDPA_poly,
                              years = 2000:2021,
                              outdir  = "mapme_data1",
                              cores = 24,
                              add_resources = TRUE)
  
  # Get GFW data
  WDPA_poly  <- get_resources(x = WDPA_poly , 
                              resources = c("gfw_treecover", "gfw_lossyear"))
  
  # Note: attr(mada_poly, "resources")$gfw_lossyear indicates the gpkg associated
  # in this case: "out_Mada/gfw_lossyear/tileindex_gfw_lossyear.gpkg"
  # So we read it to 'cut' the polygons from its borders
  footprint_treecover_tiles <- st_read(attr(WDPA_poly, "resources")$gfw_treecover) 
  footprint_lossyear_tiles <- st_read(attr(WDPA_poly, "resources")$gfw_lossyear)
} 

tmap_mode("plot")
tm_shape(footprint_treecover_tiles) +
  tm_polygons(title = "Satellite imagery", col = "papayawhip") +
  tm_shape(contour_mada) +
  tm_polygons(title = "Country") +
  tm_shape(WDPA_Mada) +
  tm_polygons(col = "green", alpha = 0.2) +
  tm_shape(footprint_treecover_tiles) +
  tm_borders(col = "red") +
  tm_add_legend(type = "fill", col = "papayawhip", 
                labels = "Source satellite imagery tiles") + 
  tm_add_legend(type = "fill", col = "grey", 
                labels = "National territory") + 
  tm_add_legend(type = "fill", col = "green", alpha = 0.3, 
                labels = "Protected areas") + 
  tm_add_legend(type = "line", col = "red", 
                labels = "No indicators if overlap") + 
  tm_layout(title = "Impossible to calculate indicators for overlapping polygons",
            title.position = c("LEFT", "TOP"))

A solution consists in segmenting all overlapping polygons along the satellite image tile borders, compute the indicators for each segment, and then aggregate the indicators for the initial polygons. Further improvement of the mapme.biodiversity package could be to implement this procedure automatically.

Code

# The following code sequence launches very heavy data processing. We want to avoid
# re-launching it if not needed. So we start checking if the output of the data 
# processing is available locally, and only if not we launch it.
if (!file.exists("data/outputs_GFC.rds")) {
  # Compute WDPA ---------------------------------------------------------------------
  
  WDPA_poly <- WDPA_poly %>%
    st_intersection(footprint_treecover_tiles) %>%
    st_make_valid() %>%
    filter(st_geometry_type(.) %in% c("POLYGON", "MULTIPOLYGON")) %>%
    # Have to re-cast back and forth because st_intersection() created 
    # multipolygons among polygons
    st_cast("MULTIPOLYGON") %>%
    st_cast("POLYGON")
  
  
  WDPA_poly <- init_portfolio(WDPA_poly,
                              years = 2000:2021,
                              outdir  = "mapme_data1",
                              cores = 24,
                              add_resources = TRUE)
  
  WDPA_poly  <- get_resources(x = WDPA_poly, 
                              resources = c("gfw_treecover", "gfw_lossyear"))
  
  WDPA_poly  <- calc_indicators(x = WDPA_poly,
                                indicators = "treecover_area", 
                                min_cover = 10, min_size = 1)
  
  WDPA_result <- WDPA_poly %>%
    select(-location) %>%
    unnest(treecover_area) %>%
    # filter(!is.na(years)) %>%
    pivot_wider(names_from = "years", values_from = "treecover", 
                names_prefix = "treecover_") %>%
    st_drop_geometry() %>%
    select(-assetid) %>%
    group_by(across(!starts_with("treecover"))) %>%
    summarise(across(starts_with("treecover"), sum, na.rm = TRUE))
  
  # Compute PAs_Vahatra -------------------------------------------------------
  
  
  PAs_Vahatra_poly <- PAs_Vahatra %>%
    st_set_crs("EPSG:4326") %>%
    st_make_valid() %>%
    st_intersection(footprint_treecover_tiles) %>%
    filter(st_geometry_type(.) %in% c("POLYGON", "MULTIPOLYGON")) %>%
    # Have to re-cast back and forth because st_intersection() created 
    # multipolygons among polygons
    st_cast("MULTIPOLYGON") %>%
    st_cast("POLYGON")
  
  
  PAs_Vahatra_poly <- init_portfolio(PAs_Vahatra_poly,
                                     years = 2000:2021,
                                     outdir  = "mapme_data1",
                                     cores = 24,
                                     add_resources = TRUE)
  
  PAs_Vahatra_poly  <- get_resources(x = PAs_Vahatra_poly, 
                                     resources = c("gfw_treecover", "gfw_lossyear"))
  
  PAs_Vahatra_poly  <- calc_indicators(x = PAs_Vahatra_poly,
                                       indicators = "treecover_area", 
                                       min_cover = 10, min_size = 1)
  
  PAs_Vahatra_result <- PAs_Vahatra_poly %>%
    select(-location) %>%
    unnest(treecover_area) %>%
    # filter(!is.na(years)) %>%
    pivot_wider(names_from = "years", values_from = "treecover", 
                names_prefix = "treecover_") %>%
    st_drop_geometry() %>%
    select(-assetid) %>%
    group_by(across(!starts_with("treecover"))) %>%
    summarise(across(starts_with("treecover"), sum, na.rm = TRUE))
  
  
  
  
  # National ------------------------------------------------------------------
  
  # Bounding box around my polygon
  bbox_mada = st_as_sf(st_as_sfc(st_bbox(contour_mada)))
  grid <- st_make_grid(x = bbox_mada,
                       cellsize = 0.3,
                       square = TRUE)
  
  mada_poly  <- contour_mada %>% 
    st_cast("POLYGON") %>%
    st_intersection(footprint_treecover_tiles) %>%
    st_intersection(grid) %>%
    st_make_valid() %>%
    filter(st_geometry_type(.) %in% c("POLYGON", "MULTIPOLYGON")) %>%
    # Have to re-cast back and forth because st_intersection() created 
    # multipolygons among polygons
    st_cast("MULTIPOLYGON") %>%
    st_cast("POLYGON")
  
  
  mada_poly <- init_portfolio(mada_poly,
                              years = 2000:2020,
                              outdir  = "mapme_data1",
                              cores = 24,
                              add_resources = TRUE)
  
  mada_poly  <- get_resources(x = mada_poly, 
                              resources = c("gfw_treecover", "gfw_lossyear"))
  
  mada_poly  <- calc_indicators(x = mada_poly,
                                indicators = "treecover_area", 
                                min_cover = 10, min_size = 1)
  
  mada_all <- mada_poly %>%
    select(-location) %>%
    unnest(treecover_area) %>%
    # filter(!is.na(years)) %>%
    pivot_wider(names_from = "years", values_from = "treecover", 
                names_prefix = "treecover_") %>%
    st_drop_geometry() %>%
    select(-assetid) %>%
    group_by(across(!starts_with("treecover"))) %>%
    summarise(across(starts_with("treecover"), sum, na.rm = TRUE))
  
  # Save the results -------------------------------------------------------
  
  save(mada_poly, WDPA_poly, PAs_Vahatra_poly, footprint_treecover_tiles,
       file = "data/outputs_GFC.rds")
  
  write_xlsx(list(Madagascar = mada_all, WDPA = WDPA_result, 
                  Vahatra = PAs_Vahatra_result),
             path = "couvert_forestier_mada_pa.xlsx")
}

4 Results

The results from this computation can be downloaded in excel by clicking this link.

As a synthesis of what precedes, the following caveats must be taken into account when interpreting this data:

• it follows the FAO (2000) definition of forest, but only imperfectly, as the minimum threshold for forest patches area had to be set to 1 hectare, instead of 0.5 hectare (a correction is forthcoming but I cannot commit on any deadline);

• the forest cover and loss estimates derivated from the WDPA data must be reported with caution, acknowledging the lack of reliability of the data source;

• the data from Vahatra seems more reliable, although it does not contain some areas reported in WDPA.

It is worth highlighting that the other information enclosed in the Vahatra data is also useful for the analysis, e.g. the representation below of protected areas date of creation, type and surface.

Code

# On ordonne les nom d'aires protégées dans l'ordre de leur séquence de création
chrono_order_PAs <- PAs_Vahatra %>%
  arrange(desc(date_creation), desc(nom)) %>%
  pull(nom)
# On transforme le champ "nom" de caractère, à une catégorisation ordonnée où
# l'ordre correspond 
PAs_Vahatra_map <- PAs_Vahatra %>%
  mutate(nom = factor(nom, levels = chrono_order_PAs),
         cat_taille = case_when(hectares > 300000 ~ 2,
                                hectares > 150000 ~ 1.5,
                                hectares >  50000 ~ 1,
                                             TRUE ~ 0.5)) %>%
  rename(`IUCN category` = cat_iucn)
# On crée un graph pour les anciennetés
graph_left <- PAs_Vahatra_map %>%
  ggplot(aes(x = date_creation, xend = ymd("2022-10-01"), y = nom, yend = nom, 
                   color = `IUCN category`)) +
  geom_segment(size = 2) +
  ggtitle("Date of creation") +
  theme(axis.ticks.y = element_blank(),
        legend.position = "none",
        axis.title = element_blank(),
        axis.text.x.bottom =  element_text(angle = 45, hjust = 1),
        axis.text.x.top = element_text(angle = -45, hjust = 1)) + 
  scale_x_date(sec.axis = dup_axis())
graph_right <- PAs_Vahatra_map %>%
  ggplot(aes(x = 0, xend = hectares/100, y = nom, yend = nom, 
                   color = `IUCN category`)) +
  geom_segment(size = 2) + 
  ggtitle("Area (km2)") +
  theme(axis.text.y = element_blank(),
        axis.title = element_blank(),
        axis.text.x.bottom =  element_text(angle = 45, hjust = 1),
        axis.text.x.top = element_text(angle = -45, hjust = 0),
        legend.position = "none") + 
  scale_x_continuous(sec.axis = dup_axis())
legende <- get_legend(graph_left  +
                        guides(color = guide_legend(nrow = 1)) +
                        theme(legend.position = "bottom"))
# On colle les deux
graphs <- plot_grid(graph_left, graph_right, rel_widths = c(2.2, 1),
          nrow = 1)
plot_grid(graphs, legende, ncol = 1,
          rel_heights = c(1,.1))

We hope that this content will be useful for the current KfW evaluation, for subsequent analysis of a similar natione and to improve the mapme.biodiversity package. For questions or comments, please contact me by email or filing an issue on the associated github repository.

Note that more comprehensive material that includes tutorials for impact evaluation and analysis is currently being prepared. The updated versions will be available here.

5 References

Chazdon, Robin L., Pedro H. S. Brancalion, Lars Laestadius, Aoife Bennett-Curry, Kathleen Buckingham, Chetan Kumar, Julian Moll-Rocek, Ima Célia Guimarães Vieira, and Sarah Jane Wilson. 2016. “When Is a Forest a Forest? Forest Concepts and Definitions in the Era of Forest and Landscape Restoration.” Ambio 45 (5): 538–50. https://doi.org/10.1007/s13280-016-0772-y.

FAO. 2000. Comparison of Forest Area and Forest Area Change Estimates Derived from FRA 1990 and FRA 2000. Forest Resources Assessment Working Paper 59. Roma. https://www.fao.org/3/ad068e/AD068E00.htm#TopOfPage.

Goodman, Steven M., Marie Jeanne Raherilalao, Sébastien Wohlhauser, Jean Clarck N. Rabenandrasana, Herivololona M. Rakotondratsimba, Fanja Andriamialisoa, and Malalarisoa Razafimpahanana. 2018. Les Aires Protégées Terrestres de Madagascar: Leur Histoire, Description Et Biote. Association Vahatra.

Görgen, Darius A., and Om Prakash Bhandari. 2022. “Mapme.biodiversity: Efficient Monitoring of Global Biodiversity Portfolios.”